function [beta nc nf]=sgd_bounded(data,lambda0,mu)
	X = data(:,1:size(data,2)-1);
	y = data(:,size(data,2));
	[n,d] = size(X);
    beta=sgd(data,0.1,0.1);
    c0=0 + 10.*rand(1,1);
    f0=-10 + 10.*rand(1,1); 
%Parameter initialization of beta
% for j=1:101
%     beta(j,:) = -10 + 20.*rand(1,d);
% 
%     if j==101
%     beta(101,:)=sgd(data,0.1,0.1);
%     end
%     for i=1:n
%         x_i = X(i,:);
%         y_i = y(i);
%         p_i = p_bounded(x_i,beta(j,:),c0,f0);
% %         yp = y_i - p_i;
%         if y_i==1
%             BLCL(i)=log(p_i);
%         else
%             BLCL(i)=log(1-p_i);
%         end
%         BLCL(i);
%     end
%     j
%     BLCLsum(j)=sum(sum(BLCL(:)));   
% end
% [a b]=find(BLCLsum==max(BLCLsum));
% figure;plot(BLCLsum);
% beta=beta(b,:);
    
	error_beta = 1;
	error_c = 1;
	error_f = 1;
 	epoch = 0;
	threshold = 10e-4;
	c = 5;
	while ((error_beta > threshold) || (error_c>threshold) || (error_f>threshold)) && epoch<40
%     for epoch=1:10
		lambda = lambda0/c^epoch;
		%lambda = lambda0/(1+lambda0*c*epoch);
		for i=1:n
                    sigmac = sigma(c0);
                    sigmaf = sigma(f0);
                    sigmanc = sigma(-c0);
                    sigmanf = sigma(-f0);
                    if sigmac==sigmaf
                        c0=c0-1;
                        f0=-f0+1;
                        continue;
                    end
                    x_i = X(i,:);
                    y_i = y(i);
                    p_i = p_bounded(x_i,beta,c0,f0);
                    yp = y_i - p_i;
            		nbeta = beta + lambda*(((yp/(p_i*(1-p_i)))*(p_i-sigmaf))*(sigmac-p_i)/(sigmac-sigmaf)*x_i-2*mu*beta);
                	nf=f0+lambda*((yp/(p_i*(1-p_i)))*sigmaf*sigmanf*(1-(1/(1+exp(-dot(x_i,beta)))))-2*mu*f0);
            		nc=c0+lambda*((yp/(p_i*(1-p_i)))*sigmac*sigmanc*(1/(1+exp(-dot(x_i,beta))))-2*mu*c0);        
            		error_beta = error_beta + sum((nbeta - beta).^2);
            		error_c = error_c + (nc-c0).^2;
            		error_f = error_f + (nf-f0).^2;
%             		error=error_beta+error_c+error_f;
            		beta = nbeta;
            		f0=nf;
            		c0=nc;
        end
        error_beta = error_beta/n;
        error_c = error_c/n;
        error_f = error_f/n;
        epoch=epoch+1
    end
     display(['Go back']);
